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A variational sine collocation method for strong-coupling problems 

Paolo AmorcQ 

Facultad de Ciencias, Universidad de Colima, 
Bernal Diaz del Castillo 340, Colima, Colima, Mexico 
(Dated: February 6, 2008) 

We have devised a variational sine collocation method (VSCM) which can be used to obtain 
accurate numerical solutions to many strong-coupling problems. Sine functions with an optimal 
grid spacing are used to solve the linear and non-linear Schrodinger equations and a lattice (j> 4 
model in (1 + 1). Our results indicate that errors decrease exponentially with the number of grid 
points and that a limited numerical effort is needed to reach high precision. 

PACS numbers: 45.10.Db,04.25.-g 
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Due to the inapplicability of perturbation theory in the 
strong coupling regime, a number of different techniques 
has been devised in the past to deal with strong-coupling 
problems. A particular attention has gone into develop- 
ing new methods, in which variational principles are used 
to improve perturbation theory, leading to results which 
are valid on a much larger domain. The Linear Delta Ex- 
pansion (LDE) Q and the Variational Perturbation The- 
ory (VPT)Q are probably the best known examples of 
such efforts: in many cases these methods allow to obtain 
series with finite (or even infinite) radius of convergence, 
in contrast with the divergent series which are usually 
obtained using perturbation theory^. 

In this letter we wish to show that the variational ideas 
which have inspired both the LDE and the VPT methods 
can be also applied to improve the performance of numer- 
ical techniques. The numerical method that we are using 
is the Sine Collocation (SCM)0, which uses sine func- 
tions to efficiently "discretize" a problem in given region. 
Sine functions are used in different areas of physics and 
mathematics ( see for example [fj and references therein) . 

A sine function is defined as 



S k (h x) = sill ( w ( x Z kh )/ h ) 
ir(x — kh)/h 

and obeys the integral representation 
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(1) 



(2) 



The reader interested in a more detailed account of the 
properties of the sine function should refer to Q ; here we 
only state the main properties which will turn useful in 
the following. 

Using eq. J2J) it is straightforward to evaluate the inte- 
grals 

/+oo 
S k (h,x)dx = h (3) 
-oo 

/ + oo 
Sk(h,x) Si(h,x) dx = h 6m- (4) 




FIG. 1: Sine functions corresponding to different values (color 
online). 



A function fix) analytic on a rectangular strip cen- 
tered on the real axis can be approximated in terms of 
sine functions as 



/(*)« J2 f(kh)S k (h,x) 



(5) 



Using eq. © together with eq. © one obtains 

/+oo oo 
/(*) dx w h Yl f( k > h )- ( fi ) 

K — — OQ 

An expression for the error in eq. (JSJ) has been obtained 
by Stenger 0, showing that it decays exponentially as 
the spacing h is reduced. 

As the reader can appreciate in FigQ for a fixed h a 
given sine function selects a point on the real line, cor- 
responding to its maximum value, and vanishes in the 
points of intersection with the other sine functions. This 
property, which allows to obtain a proper "discretiza- 
tion" of a problem in the continuum, is at the basis of 
the SCM. 

We now describe the SCM in detail by considering the 
stationary Schrodinger equation 



V(x^hlx) = Eib(r) . 



(7) 
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affects the precision of the numerical results, no proce- 
dure to determine h is discussed in 0. Using a different 
method, the author and collaborators Q have solved nu- 
merically the Schrodinger equation for the anharmonic 
oscillator using an arbitrary basis of Gauss-Hermitc func- 
tions, depending upon a scale factor. In that paper it was 
proved that the arbitrary scale factor can be chosen op- 
timally by applying the principle of minimal sensitivity 
(PMS)H to the sub-trace of the hamiltonian matrix. 

Using the same procedure we regard h as a variational 
parameter and consider the trace 



FIG. 2: log 10 \Eo — Eo xact \ as a function of the spacing h using 



10 for the harmonic oscillator V{x) = x /2. 



Tr l H ] = ( 2k ma* + 1) + J2 V ^ ' C 11 ) 



k=-k n 



The matrix elements Hki of the hamiltonian evaluated 
in the set of sine functions are given by 



where 2k max + 1 is the number of sine functions (grid 
points) used in the evaluation. 

The solution to the PMS equation 



H 



hi 



-\ c u+^i v ( kh ) 



(8) 



Tk Tr[H] 



(12) 



Notice that the kinetic term has been obtained by using 
the property 



d 2 °° 
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where 



(2) 
c lk 
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C ik S k (h,x) , 



if k = l 



(9) 



(10) 



while the potential matrix has been approximated by the 
diagonal matrix of the potential evaluated over the grid. 

Once h is specified the diagonalization of H^i allows 
to obtain numerical approximations to the energies and 
wave functions of the problem: this strategy was used 
in to solve the Schrodinger equation corresponding to 
different potentials. Although the choice of h strongly 



provides the optimal spacing [12]]. Once that h has been 
determined using the PMS our method allows to obtain 
quite rapidly the numerical approximations [13j. 

In fig. [5] we display the log 10 \Eq — EQ Xact \ as a function 
of the spacing h using 21 grid points (k max = 10) for the 
harmonic oscillator V{x) = x 2 /2. The PMS condition, 
cq. H12(l . yields hp^s = 0.547, which is remarkably close 
to the minimum of the curve. 

The anharmonic oscillator 



H = -— +x z +g x q 



provides a more demanding test of our method. We 
have obtained the ground state energy corresponding to 
g = 2000 using a grid of 101 points (kmax = 50) and 
compared it with the precise results of |9|, finding that 
the first 42 digits are correct (underlined): 



E a = 13.38844170100806193900617690280728652296099 . 



The result is also seen to converge exponentially to the 
exact answer as a function of the number of grid points. 
We now consider the Gross-Pitaevskii (GP) equation 



1 d*_ 

2 dx 2 



+ V ext {x) +4 Tr a \ip{x){' 



ip(x) = E ip(x) 



which is relevant in the study of Bose-Einstein conden- 

Qtttinn TA^f-r^ i« tVin OYtornwl ennfini-nev nntnntinl flnrl 



E is the energy of the condensate. The wave function 
ip{x) is normalized to yield the number of particles in 
the condensate, i.e. f dx |?/>(a;)| 2 = N. 

We have applied our method to the GP equation by 
first solving the corresponding linear equation, and by 
then implementing a self-consistent procedure in which 
the density term is evaluated taking the wave function 

rnlrnlfltprl at tVin nrmnmiQ Qtnn MnH ttinn it iq iiqpH tn KiiilH 
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FIG. 3: Probability density for the ground state of the GP 
equation using V ex t(x) = x 1 11 (color online). 



tion corresponds to the usual harmonic oscillator wave 
function. In fig. 0] we have plotted log 10 \E — Eo xact \ as 
a function of the number of iterations for different grid 
sizes (11, 21 and 31 respectively): the datas initially dis- 
play an exponential decay - independent of the grid size 
- which is then followed by a plateau. The plateau sig- 
nals that the maximal precision has been achieved for a 
given grid size. 

As a last example of application of our method we 
consider a lattice <fi 4 in 1 + 1 dimensions. This model has 
been studied by Nishiyama in [lfl | and corresponds to 
the hamiltonian 
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(13) 
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The fields obey the canonical commutation relations 
[4>i,-Kj] = iSij and [4>i,4>j] = [-Ki^j] = 0. Following 01 
we perform a rescaling of the fields — > g _1 / 6 and 
and obtain 



H = g 



1/3 



E 



2 



l+ i) 2 + 



A2\ 1 



iterations 



FIG. 4: log 10 \Eq — Ef ) xact \ as a function of the number of 
iterations for different grid sizes (color online). 



an effective potential V e ff(x) = V ex t(x) + 4 it a \tp(x)\ 2 . 
In this potential the resulting Schrodingcr equation is 
solved again and the procedure is iterated until self- 
consistency is reached. 

In fig. |3 we have plotted the wave function in an har- 
monic trap, obtained after 0, 20 and 30 iterations of our 
method, assuming Ana = 1 and \ip(x)\ 2 dx = 2. A 
grid of 21 points has been used. The 0th order wave func- 



where A = g 2 / 3 . We express the ground state energy as 

TP. — a 1 / 3 f 

Nishiyama has numerically solved this model using a 
Linked Cluster Expansion (LCE) and the Density Ma- 
trix Renormalization Group (DMRG) llj : using the LCE 
he has obtained a perturbation series in A, up to or- 
der 11, whose convergence has then been improved using 
Aitken's 6 2 process. The comparison with the DMRG 
results shows that LCE is valid up to A « 2. 

We wish to show that the same problem can be 
solved using our method. We have proceeded as fol- 
lows: first we have solved the Schrodinger equation 
for the anharmonic oscillator, corresponding to setting 
A = 0, and we have obtained the wave function ^(4>) ~ 
Sr=- X fe a r S r (h, </>); we have then used ^f(4>) to evalu- 
ate the matrix element (<I'(</)i+i)|i/| , I , (</);+i)), obtaining 
the effective potential felt by the ith site: 



V{ct>) = 4 + A 



2<j) £ a l rh + £ a 2 r {rh) 2 



(14) 



=-fe. 



The Schrodingcr equation is then solved again and the 
coefficients a r are recalculated. The procedure is re- 
peated until self-consistency is reached. 

In fig. we have calculated the scaled ground state en- 

prefv t= _ fi« a funrtinn nf \ llfiincr a crriH nf 41 crriH r>nintQ 



(kmax = 20). The solid curve corresponds to the pertur- 
bative expansion of eq.(6) of 0], to order A 11 , obtained 
with the LCE. Our result compares quite favorably also 
with the results obtained with the DMRG - sec Fig. 2 of 
Hoi and have been obtained in few minutes of running 
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FIG. 5: Ground state energy for the lattice <^> 4 model as a 
function of A = a~ 2//3 (color online). 



grid spacings are obtained using the PMS and allow to 
achieve optimal results with limited numerical effort; the 
diagonalization of the hamiltonian matrix provides ap- 
proximations for the energies and wave functions of part 
of the spectrum and can used to study the time evolution 
of a wave packet (similarly to what done in Q); finally, 
our method can also be applied without modifications to 
models with non-polynomial interactions. 

The last example shows that VSCM can be a useful 
tool in the numerical solution of lattice quantum models 
and it could possibly provide an alternative to numerical 
methods already present on the "market" . Future work 
in this direction is expected. 

P.A. acknowledges support of Conacyt grant no. C01- 
40633/A-l. 




FIG. 6: Effective potential V(tf>) for values of A close to the 
discontinuity. The curves differ by AA = 0.1. 



time on a Linux desktop running Mathcmatica. 

Notice that in [ljj it was speculated the presence of 
a singularity around A = — 2, which is possibly related 
to the onset of a phase transition. Our simulation shows 
the presence of a discontinuity located at A w —1.75. As 
it can be seen from fig. such discontinuity is due to a 
sudden jump in the effective potential. 

We wish to conclude this letter stressing few points, 
which we believe to be important: the VSCM provides 
errors which decrease exponentially with the grid size; the 
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